A data-driven approach to evaluating landslide hazards
Elsa Culler
Data
This study uses landslide, fire, and precipitation data to study the antecedent precipitation of post-wildfire landslides. The impact of wildfire on rainfall-triggered landslide susceptibility in a variety of global regions is evaluated through a proxy of normalized precipitation.
Data Filepaths
Use this section to modify file paths if running the notebook on a different computer.
# Landslide data from the NASA Global Landslide Catalog
slide.path <- '../regions_data/glc/GLC20180821.csv'
# Fire data from the MODIS Burned Area dataset
modis.paths <- list.files(path='../regions_data/processed',
pattern='slide_modisburndate_[ns][ew].csv',
full.names = T)
# Precipitation data from CHIRPS
precip.chirps.ww.path <- file.path("../regions_data/processed",
"glc_precip_global.csv")
precip.chirps.ww.monthly.path <- file.path("../regions_data/", "processed",
"glc_precip_global_monthly.csv")
# Additional preciptiation data from Daymet for a limited region
precip.daymet.ca.path <- file.path("../regions_data/processed",
"glc_precip_daymet.csv")Data inclusion
Include landslides with the following parameters: * Rainfall triggered
rain.triggers <- c('rain', 'downpour', 'flooding', 'continuous_rain')- Occurred within the time frame of the fire and precipitation datasets
start.date <- lubridate::ymd('2000-01-01')
end.date <- lubridate::ymd('2019-01-01')- Occurred within the spatial extend of the fire and precipitation datasets
min.latitude <- -50
max.latitude <- 50- Location accuracy on par with precipitation resolution or more accurate
min.location.accuracy <- '10km'- MODIS Burned Area data should include an additional 3 years
modis.start.date <- lubridate::ymd('20040101')
modis.end.date <- lubridate::ymd('20190101')- Wilcox tests between burned and unburned locations have the following parameters:
wilcox.alternative <- 'less' # n - burned, m - unburned
wilcox.alpha <- 0.05Study definitions
- Study location parameters are applied using a set of functions to identify if data are included in the study.
# Identify if a location is in the global limits
is.in.global.study <- function (x, y) { between(y, min.latitude, max.latitude) }
# Identify if a location is in the CA/NV pilot study
get.data('state.boundaries', tigris::states)
boundaries <- subset(state.boundaries, STUSPS %in% c('CA', 'NV'))
boundaries.df <- fortify(boundaries)
is.in.ca.study <- function (x, y) {
bbox <- sf::st_bbox(subset(state.boundaries, STUSPS == 'CA'))
between(x, bbox$xmin, bbox$xmax) &
between(y, bbox$ymin, bbox$ymax)
}Seasons are split in the Northern Hemisphere by:
- Winter
- December-January-February
- Spring
- March-April-May
- Summer
- June-July-August
- Fall
- September-October-November
In the southern hemisphere, the split is shifted instead:
- Winter
- June-July-August
- Spring
- September-October-November
- Summer
- December-January-February
- Fall
- March-April-May
A set of functions converts months and days-of-the-year to seasons.
season.breaks <- c('spring'=59, 'summer'=151, 'fall'=243, 'winter'=334)
yday.to.season <- function(yday, latitude) {
if (latitude > 0) {
case_when(
yday > season.breaks['winter'] ~ 'Winter',
yday > season.breaks['fall'] ~ 'Fall',
yday > season.breaks['summer'] ~ 'Summer',
yday > season.breaks['spring'] ~ 'Spring',
TRUE ~ 'Winter')
} else {
case_when(
yday > season.breaks['winter'] ~ 'Summer',
yday > season.breaks['fall'] ~ 'Spring',
yday > season.breaks['summer'] ~ 'Winter',
yday > season.breaks['spring'] ~ 'Fall',
TRUE ~ 'Summer')
}
}
month.to.season <- function (month, north) {
if (north) {
case_when(
month >= 12 ~ 'Winter',
month >= 9 ~ 'Fall',
month >= 6 ~ 'Summer',
month >= 3 ~ 'Spring',
TRUE ~ 'Winter')
} else {
case_when(
month >= 12 ~ 'Summer',
month >= 9 ~ 'Spring',
month >= 6 ~ 'Winter',
month >= 3 ~ 'Fall',
TRUE ~ 'Summer')
}
}- Consistent plot colors and themes
# Region colors
region.plot.order <- c('All landslides',
'California',
'Intermountain West',
'Pacific Northwest',
'Central America',
'Himalayas',
'Southeast Asia')
scale_region_values <- c('All landslides' = 'grey50',
setNames(RColorBrewer::brewer.pal(6, 'Dark2'),
region.plot.order[2:7]))
scale_region_labels <- c('All landslides' = 'None',
setNames(region.plot.order[2:7],
region.plot.order[2:7]))
scale_color_region <- scale_color_manual('Region',
values = scale_region_values,
labels = scale_region_labels)
scale_fill_region <- scale_fill_manual('Region',
values = scale_region_values,
labels = scale_region_labels)
scale_color_precip <- scale_color_distiller('Precipitation percentile',
palette='BrBG', direction = 1,
guide = guide_colorbar(barwidth = unit(0.5, 'npc')))
scale_fill_precip <- scale_fill_distiller('Precipitation percentile',
palette='BrBG', direction = 1,
guide = guide_colorbar(barwidth = unit(0.5, 'npc')))
# Fire colors
burned.labels <- c('burned'='Burned', 'unburned'='Unburned')
burned.colors <- c('burned' = '#d95f02', unburned = '#7570b3')
scale_color_burned <- scale_color_manual('',
values = burned.colors,
labels = burned.labels)
scale_fill_burned <- scale_fill_manual('',
values = burned.colors,
labels = burned.labels)
# Season
scale_x_season <- scale_x_continuous(
breaks = season.breaks,
labels = c('Spring', 'Summer', 'Fall', 'Winter'))
col.ocean <- 'lightskyblue1'theme_map <- list(
theme_bw(),
coord_quickmap(),
labs(x='', y='')
)
theme_global_map <- list(theme_map,
borders(fill = 'grey90', size = 0.1),
theme(panel.background = element_rect(fill=col.ocean),
panel.grid = element_blank(),
axis.text = element_blank(),
axis.ticks = element_blank(),
plot.margin = unit(c(0,0,0,0), "null")),
coord_quickmap(ylim = c(-60, 75), xlim = c(-180, 180))
)
theme_panel <- theme(panel.background = element_rect(fill=col.ocean),
panel.grid = element_blank(),
axis.text = element_blank(),
axis.ticks = element_blank(),
plot.margin = unit(c(1,0,-1,0), "lines"),
legend.position = 'null')
theme_panel_legend <- theme(
legend.box.margin = margin(-1, 0, 1, 0),
legend.box.just = c(0, 0),
#legend.background = element_rect(color='black'), #debug
legend.position = 'bottom'
)
capitalize <- function(string) {
substr(string, 1, 1) <- toupper(substr(string, 1, 1))
string
}Landslide data from the NASA Global Landslide Catalog
load.slides <- function (slide.path,
start.date, end.date,
is.in.global.study,
rain.triggers,
min.location.accuracy) {
# Define factor levels
ls.sizes <- c("small", "medium", "large", "very_large", "catastrophic")
ls.accuracy <- c("exact", "1km", "5km", "10km",
"25km", "50km", "100km", "250km")
# Read data file and remove duplicates
slide.raw <- readr::read_csv(slide.path,
guess_max = 100000,
na = c("", "unknown", "--"),
col_types = cols(
landslide_size = readr::col_factor(ls.sizes, ordered=TRUE),
location_accuracy = readr::col_factor(ls.accuracy, ordered=TRUE),
landslide_trigger = readr::col_factor(NULL),
landslide_category = readr::col_factor(NULL),
landslide_setting = readr::col_factor(NULL),
country_code = readr::col_factor(NULL)
)
) %>%
dplyr::mutate(event_date = as_date(ymd_hm(event_date)),
submitted_date = as_date(ymd_hm(submitted_date)),
last_edited_date = as_date(ymd_hm(last_edited_date)))
slide.unique <- slide.raw %>%
dplyr::group_by(event_date, latitude, longitude) %>%
dplyr::mutate(i = row_number()) %>%
dplyr::filter(i == 1) %>%
dplyr::select(-i) %>%
dplyr::ungroup()
# Filter landslides to meet study parameters
slide.unique %>%
dplyr::filter(event_date >= start.date,
event_date <= end.date,
is.in.global.study(longitude, latitude),
landslide_trigger %in% rain.triggers,
location_accuracy <= min.location.accuracy) %>%
droplevels() %>%
dplyr::mutate(event_yday = yday(event_date),
debris_flow = landslide_category %in% c('mudslide', 'debris flow'),
season = factor(map2_chr(event_yday, latitude,
~yday.to.season(.x, .y)),
levels = c('Spring', 'Summer', 'Fall', 'Winter'),
ordered = T))
}
get.data('slide.all.precip', overwrite = F,
load.slides,
slide.path, start.date, end.date,
is.in.global.study, rain.triggers,
min.location.accuracy)
slide.all.precip %>%
dplyr::select(landslide_trigger, location_accuracy, latitude, event_date) %>%
summary()## landslide_trigger location_accuracy latitude
## rain :1823 exact: 466 Min. :-46.77
## continuous_rain: 588 1km :1707 1st Qu.: 10.49
## downpour :3254 5km :2437 Median : 29.91
## flooding : 48 10km :1103 Mean : 23.52
## 3rd Qu.: 39.11
## Max. : 49.98
## event_date
## Min. :2000-10-13
## 1st Qu.:2010-11-30
## Median :2013-06-24
## Mean :2013-05-17
## 3rd Qu.:2016-03-11
## Max. :2018-08-07Precipitation data from CHIRPS and Daymet
Global CHIRPS data
load.precip.all <- function (precip.chirps.ww.path) {
precip.chirps.ww.raw <- read_csv(precip.chirps.ww.path,
col_types = cols(X1 = col_datetime(format = "%Y-%m-%d")))
precip.chirps.ww.raw %>%
dplyr::select(date=X1, OBJECTID, pctl=precip_pctl_7) %>%
dplyr::filter(OBJECTID %in% slide$OBJECTID) %>%
dplyr::mutate(window=7) %>%
dplyr::filter(!is.na(pctl)) # Remove NA values from beginning of data range
}
get.data('precip.chirps.ww.all', precip.chirps.ww.path)
precip.chirps.ww.allDetailed global CHIRPS data close to the landslide
date.filter.precip <- function (raw, slide=slide.all.precip, pre=90, post=30) {
raw %>%
dplyr::right_join(slide %>% dplyr::select(OBJECTID, event_date),
by='OBJECTID') %>%
dplyr::select(OBJECTID, dplyr::starts_with('precip_'), event_date, date) %>%
dplyr::mutate(date = lubridate::as_date(date),
days_before = as.integer(event_date - date)) %>%
dplyr::filter(days_before < pre, days_before > -post)
}
tidy.precip <- function (precip) {
precip %>%
dplyr::select(OBJECTID, starts_with('precip_'), event_date, date) %>%
gather(key='variable', value='amount', starts_with('precip_')) %>%
separate(variable, into=c(NA, 'variable', 'window')) %>%
dplyr::filter(variable=='pctl') %>%
dplyr::select(-variable, pctl=amount) %>%
dplyr::mutate(days_before = as.integer(event_date - date),
yday = yday(event_date - days(days_before)),
window = as.numeric(window)) %>%
dplyr::select(-event_date)
}
load.event.precip <- function (precip.chirps.ww.path) {
precip.chirps.ww.raw <- read_csv(precip.chirps.ww.path,
col_types = readr::cols(X1 = col_datetime(format = "%Y-%m-%d")))
precip.chirps.ww.raw.filt <- precip.chirps.ww.raw %>%
dplyr::rename(date = X1) %>%
dplyr::filter(year(date) > 2003)
precip.chirps.ww.raw.filt2 <- precip.chirps.ww.raw.filt %>%
date.filter.precip
precip.chirps.ww.raw.filt2 %>%
tidy.precip
}
get.data('precip.chirps.ww.event', precip.chirps.ww.path)
precip.chirps.ww.eventImport Monthly CHIRPS data
load.precip.monthly <- function (monthly.path, slide) {
read_csv(precip.chirps.ww.monthly.path,
col_types = cols(X1 = col_date(format = '%Y-%m-%d'))) %>%
dplyr::transmute(OBJECTID,
year = year(X1), month = month(X1),
precip.mm=precip) %>%
dplyr::inner_join(slide %>% select(OBJECTID, latitude), by='OBJECTID') %>%
dplyr::mutate(north = latitude > 0) %>%
dplyr::group_by(north) %>%
dplyr::mutate(season = factor(month.to.season(month, north),
levels = c('Spring', 'Summer', 'Fall', 'Winter'),
ordered = T)) %>%
ungroup()
}
get.data('precip.chirps.ww.monthly', overwrite = F,
load.precip.monthly,
precip.chirps.ww.monthly.path,
slide.all.precip)
precip.chirps.ww.seasonal.median <- precip.chirps.ww.monthly %>%
group_by(OBJECTID, year, season) %>%
summarize(precip.mm = sum(precip.mm)) %>%
ungroup() %>%
summarize(median = median(precip.mm)) %>%
pull(median)
precip.chirps.ww.monthlyFilter out landslides with zero precipitation
precip.chirps.ww.event %>%
# Filter out sites with no precipitation data at all
filter(days_before==-1, window==7, !is.na(pctl)) %>%
# Filter out sites already excluded by landslide criteria
inner_join(slide.all.precip, by = 'OBJECTID') %>%
# Count numbers of zero and non-zero precipitation
group_by(precip = ifelse(pctl > 0, 'non.zero.precip', 'zero.precip')) %>%
tally() %>%
spread(key = precip, value = n) %>%
mutate(total = non.zero.precip + zero.precip,
percent.zero.precip = zero.precip / total)filter.slide.precip <- function (slide.all.precip, precip.event) {
slide.all.precip %>%
dplyr::inner_join(precip.event %>%
filter(days_before==-1, window==7, !is.na(pctl)) %>%
dplyr::filter(days_before==-1, window==7, pctl > 0)) %>%
dplyr::select(OBJECTID, latitude, longitude, event_date, event_yday, season,
location_accuracy,
landslide_category, landslide_trigger, landslide_size,
landslide_setting, debris_flow)
}
get.data('slide', overwrite = F,
filter.slide.precip,
slide.all.precip,
precip.chirps.ww.event)
slideFilter California Landslides
filter.slide.ca <- function (slide.all.precip) {
slide.all.precip %>%
filter(is.in.ca.study(longitude, latitude)) %>%
mutate(region = as.factor(case_when(
latitude < 35.0 ~ 'South',
longitude < -121.25 ~ 'North',
TRUE ~ 'East')))
}
get.data('slide.ca', overwrite = F,
filter.slide.ca,
slide.all.precip)
slide %>% summary## OBJECTID latitude longitude event_date
## Min. :388907 Min. :-46.77 Min. :-179.865 Min. :2007-01-05
## 1st Qu.:391763 1st Qu.: 10.33 1st Qu.: -96.159 1st Qu.:2010-11-30
## Median :394578 Median : 28.55 Median : -8.624 Median :2013-06-16
## Mean :394614 Mean : 23.09 Mean : -2.247 Mean :2013-05-15
## 3rd Qu.:397488 3rd Qu.: 38.95 3rd Qu.: 92.212 3rd Qu.:2016-03-10
## Max. :400281 Max. : 49.98 Max. : 179.402 Max. :2018-08-07
##
## event_yday season location_accuracy landslide_category
## Min. : 1.0 Spring:1310 exact: 423 landslide :3492
## 1st Qu.: 91.0 Summer:1798 1km :1573 mudslide :1326
## Median :172.0 Fall :1104 5km :2259 rock_fall : 226
## Mean :171.6 Winter:1101 10km :1058 complex : 103
## 3rd Qu.:244.0 debris_flow: 98
## Max. :366.0 (Other) : 64
## NA's : 4
## landslide_trigger landslide_size landslide_setting
## rain :1656 small :1735 above_road :1257
## continuous_rain: 561 medium :3160 natural_slope: 265
## downpour :3054 large : 326 urban : 178
## flooding : 42 very_large : 36 below_road : 121
## catastrophic: 3 above_river : 65
## NA's : 53 (Other) : 187
## NA's :3240
## debris_flow
## Mode :logical
## FALSE:3987
## TRUE :1326
##
##
##
## Filter California CHIRPS data
filter.precip.ca <- function (precip.chirps.ww.event, slide.ca) {
precip.chirps.ca <- precip.chirps.ww.event%>%
dplyr::filter(OBJECTID %in% slide.ca$OBJECTID)
}
get.data('precip.chirps.ca.event', filter.precip.ca, overwrite = F,
precip.chirps.ww.event, slide.ca)
get.data('precip.chirps.ca.all', filter.precip.ca, overwrite = F,
precip.chirps.ww.all, slide.ca)
precip.chirps.ca.event %>%
# Filter out sites with no precipitation data at all
filter(days_before==-1, window==7, !is.na(pctl)) %>%
# Filter out sites already excluded by landslide criteria
inner_join(slide.all.precip, by = 'OBJECTID') %>%
# Count numbers of zero and non-zero precipitation
group_by(precip = ifelse(pctl > 0, 'non.zero.precip', 'zero.precip')) %>%
tally() %>%
spread(key = precip, value = n) %>%
mutate(total = non.zero.precip + zero.precip,
percent.zero.precip = zero.precip / total * 100)Import Daymet data for California
load.precip.daymet.ca <- function (precip.daymet.ca.path) {
read_csv(precip.daymet.ca.path,
col_types = cols(X1 = col_datetime(format = "%Y-%m-%d %H:%M:%S"))) %>%
rename(date=X1) %>%
date.filter.precip %>%
tidy.precip
}
get.data('precip.daymet.ca.all', overwrite = F,
load.precip.daymet.ca,
precip.daymet.ca.path)
precip.daymet.ca.allFire data computed from MODIS Burned Area
load.modis <- function (modis.paths, burn.date.min, burn.date.max) {
modis.date.raw <- dplyr::bind_rows(map(modis.paths, read_csv))
burn.record.length <- time_length(burn.date.max - burn.date.min, 'years')
modis.date.raw %>%
dplyr::right_join(slide %>% dplyr::select(OBJECTID, event_date, latitude),
by='OBJECTID') %>%
# Convert MODIS Burned Area value to date
dplyr::mutate(burn_date = make_date(year=floor(month / 1000)) + burn_date) %>%
dplyr::select(-month) %>%
# Separate out indivicual fires
dplyr::group_by(OBJECTID) %>%
dplyr::mutate(gap = as.integer(burn_date - lag(burn_date, default=min(burn_date))),
new.fire = gap > 30,
fire.id = cumsum(new.fire)) %>%
drop_na(fire.id) %>%
# Calculate burn fraction, fire ending date, and season for each fire
dplyr::group_by(OBJECTID, event_date, fire.id) %>%
summarize(fraction = sum(fraction),
burn_date = max(burn_date),
burn_season = yday.to.season(yday(burn_date), latitude)) %>%
# Calculate time relative to landslide
dplyr::mutate(burn_days_before = event_date - burn_date,
burn_years_before = time_length(burn_days_before, 'years')) %>%
# Compute fire information
dplyr::group_by(OBJECTID) %>%
dplyr::mutate(pre.slide = burn_date==max(burn_date[burn_days_before>=0]),
n = n(),
period = burn.record.length / n,
start_days_before = min(event_date) - burn.date.min,
start_years_before = time_length(start_days_before, 'years'),
end_days_before = max(event_date) - burn.date.max,
end_years_before = time_length(end_days_before, 'years')) %>%
dplyr::select(-event_date) %>%
dplyr::ungroup()
}
get.data('modis.date', overwrite = F,
load.modis,
modis.paths, modis.start.date, modis.end.date)
modis.dateclassify.burned <- function (modis.date, slide) {
modis.date %>%
dplyr::filter(burn_years_before < 3, pre.slide) %>%
dplyr::right_join(slide, by='OBJECTID') %>%
dplyr::group_by(OBJECTID) %>%
summarize(burn.cumulative = sum(fraction)) %>%
dplyr::mutate(burn.cumulative = ifelse(is.na(burn.cumulative),
0,
burn.cumulative),
burned = ifelse(burn.cumulative > 0, 'burned', 'unburned')) %>%
dplyr::select(OBJECTID, burned, burn.cumulative) %>%
dplyr::ungroup()
}
get.data('modis', overwrite = F,
classify.burned,
modis.date,
slide)
modis %>%
group_by(burned) %>%
tally() %>%
mutate(percent = n / nrow(slide) * 100)Define regions using AGNES clustering
cluster.global <- function (slide, agclusts, k=30) {
# Select clusters with more than 100 landslides, and label
slide %>%
dplyr::select(OBJECTID, latitude, longitude) %>%
dplyr::mutate(region = agclusts %>% cutree(k)) %>%
dplyr::group_by(region) %>%
mutate(region = as_factor(region)) %>%
dplyr::filter(n() > 100) %>%
dplyr::ungroup() %>%
dplyr::mutate(region = fct_recode(as_factor(region),
'Central Europe' = '1',
'Malaysia' = '2',
'Western US/Canada' = '3',
'East Brazil' = '4',
'Southeast Asia' = '5',
'South India' = '6',
'Eastern US' = '8',
'Himalayas' = '9',
'Southeast Asia' = '10',
'Caribbean/Venezuela' = '11',
'Central America' = '14'))
}
# Cluster landslides based on location using AGNES algorithm
get.data('agclusts',
cluster::agnes,
slide %>%
dplyr::select(latitude, longitude))
get.data('clusters.global', overwrite = F,
cluster.global,
slide.all.precip, agclusts)
clusters.global %>%
ggplot(aes(longitude, latitude)) +
theme_global_map +
geom_point(aes(color = region), size=1) +
scale_color_brewer('Region', palette='Paired') +
guides(color=guide_legend(override.aes = list(size=5)))split.regions.us <- function (agclusts, slide, k=6) {
# Crop cluster tree, select clusters with more than 100 landslides, and label
slide %>%
dplyr::select(OBJECTID, latitude, longitude) %>%
dplyr::mutate(region = agclusts %>% cutree(k)) %>%
mutate(region = as_factor(region)) %>%
dplyr::mutate(region = fct_recode(as_factor(region),
'Intermountain West' = '3',
'Intermountain West' = '4',
'Intermountain West' = '5',
'Pacific Northwest' = '1',
'California' = '2')) %>%
dplyr::group_by(region) %>%
dplyr::filter(n() > 100) %>%
dplyr::ungroup() %>%
droplevels()
}
get.data('agclusts.us', overwrite = F,
function (coord) cluster::agnes(coord, stand=T),
clusters.global %>%
dplyr::filter(region=='Western US/Canada') %>%
dplyr::select(latitude, longitude))
get.data('clusters.us', overwrite = F,
split.regions.us,
agclusts.us,
clusters.global %>%
dplyr::filter(region=='Western US/Canada'))
clusters.us %>%
ggplot(aes(longitude, latitude)) +
geom_point(aes(color = region), size=1) +
scale_color_brewer(palette='Set2') +
guides(color=guide_legend(override.aes = list(size=5))) +
theme_mapmerge.regions <- function (clusters.global, clusters.us) {
# Merge US and global regions
clusters.global %>%
left_join(clusters.us %>% dplyr::select(OBJECTID, region),
by='OBJECTID') %>%
dplyr::mutate(region = as_factor(ifelse(!is.na(region.y),
as.character(region.y),
as.character(region.x)))) %>%
dplyr::select(-starts_with('region.')) %>%
filter(region != 'Western US/Canada') %>%
droplevels()
}
get.data('regions.all', overwrite = F,
merge.regions,
clusters.global, clusters.us)get.regions.info <- function (regions.all, slide) {
regions.all %>%
bind_rows(slide %>% dplyr::mutate(region = 'All landslides')) %>%
dplyr::inner_join(modis, by = 'OBJECTID') %>%
group_by(region) %>%
dplyr::group_by(region, burned) %>%
dplyr::summarize(n = n()) %>%
tidyr::spread(key = burned, value = n) %>%
dplyr::mutate(n = burned + unburned,
percent.burned = burned / n * 100) %>%
ungroup
}
get.data('regions.info', overwrite = F,
get.regions.info,
regions.all, slide)
regions.infofilter.and.merge.regions <- function (regions.all, regions.info, slide) {
# Combine nearby regions
regions.all %>%
filter(OBJECTID %in% slide$OBJECTID) %>%
bind_rows(slide %>% dplyr::mutate(region = 'All landslides')) %>%
dplyr::left_join(regions.info) %>%
mutate(region = fct_recode(region,
'Central America' = 'Caribbean/Venezuela')) %>%
dplyr::filter(region != 'East Brazil',
percent.burned > 4) %>%
droplevels() %>%
dplyr::mutate(region = factor(region, region.plot.order, ordered = T)) %>%
dplyr::select(OBJECTID, region)
}
get.data('slide.regions', overwrite = F,
filter.and.merge.regions,
regions.all,
regions.info,
slide)
regions.tally <- slide.regions %>%
group_by(region) %>%
tally
regions.tallyExploration
Map
A panelled maps showing landslide regions (a and b), whether or not each landslide took place on a burned site (c and d), and the precipitation percentile of the landslide-triggering storm (e and f). Inset plots show, for each region, the distribution of burned fractions of the landslide sites (a and b), location accuracy of landslides split by fire history (c and d), and the seasonal total precipitation climatology (e and f). ### Regions map
gg.region <- slide.regions %>%
left_join(slide, by='OBJECTID') %>%
mutate(region = fct_relevel(region, 'All landslides', after = Inf)) %>%
arrange(desc(region)) %>%
ggplot(aes(longitude, latitude)) +
theme_global_map +
geom_point(aes(color = region), size=0.02) +
labs(color = '') +
scale_color_region +
guides(color = guide_legend(override.aes = list(size = 2))) +
theme_panel_legend
gg.regionBurned sites map
gg.burned <- slide %>%
dplyr::left_join(modis) %>%
dplyr::arrange(desc(burned)) %>%
ggplot() +
theme_global_map +
#geom_point(data = slide.regions %>%
# dplyr::filter(region !='All landslides') %>%
# dplyr::left_join(slide, by='OBJECTID'),
# aes(x = longitude, y = latitude, fill = region),
# size = 2, alpha = .01, shape = 21, stroke = NA) +
geom_point(aes(x=longitude, y=latitude, color=burned),
size=.02) +
scale_color_burned +
#scale_fill_region +
guides(color = guide_legend(override.aes = list(size = 2),
title = NULL),
fill = guide_none()) +
theme_panel_legend
gg.burnedPrecipitation map
gg.precip <- slide %>%
left_join(precip.chirps.ww.event %>%
filter(days_before==-1, window==7),
by='OBJECTID') %>%
filter(pctl > 0) %>%
arrange(desc(pctl)) %>%
ggplot() +
theme_global_map +
geom_point(aes(x=longitude, y=latitude, color=pctl),
size=.02) +
scale_color_precip +
theme_panel_legend
gg.precipBurned fraction
min.burn.cumulative <- modis %>%
filter(burn.cumulative > 0) %>%
pull(burn.cumulative) %>%
min
max.burn.cumulative <- modis %>%
pull(burn.cumulative) %>%
max
plot.burned.fraction <- function (reg, slide.regions, modis) {
slide.regions %>%
filter(region==reg) %>%
left_join(modis, by = 'OBJECTID') %>%
filter(burn.cumulative > 0) %>%
ggplot() +
geom_density(aes(x = burn.cumulative), fill = 'black') +
scale_x_log10() +
labs(x = 'Burned fraction', y = '') +
scale_color_region +
scale_fill_region +
coord_cartesian(expand=F, xlim = c(min.burn.cumulative, max.burn.cumulative)) +
theme_bw() +
theme(legend.position = 'none',
axis.text.x = element_text(angle = 90, vjust = 0.5, size = 4),
axis.title = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
plot.background = element_blank(),
panel.grid = element_blank(),
plot.title = element_text(size = 5, margin=margin(0,0,0,0))) +
ggtitle(str_c(reg, '\nn=', pull(filter(regions.tally, region==reg))))
}
gg.burned.fraction <- region.plot.order %>%
set_names %>%
map(plot.burned.fraction, slide.regions, modis)
cowplot::plot_grid(plotlist = gg.burned.fraction, nrow = 1)Location accuracy by region/fire history
plot.accuracy <- function (reg, slide.regions, slide, modis) {
slide.regions %>%
left_join(slide, by='OBJECTID') %>%
left_join(modis, by='OBJECTID') %>%
filter(region==reg) %>%
mutate(location_accuracy = fct_recode(location_accuracy, '0km' = 'exact')) %>%
drop_na(location_accuracy) %>%
ggplot() +
geom_bar(aes(x = location_accuracy, fill=burned,
group = interaction(region, burned),
y = ..prop..),
position=position_dodge2(preserve = 'single')) +
scale_fill_burned +
coord_cartesian(expand=F, ylim = c(0, 0.6)) +
theme_bw() +
labs(y = '', x = '', fill = 'Location\naccuracy') +
theme(legend.position = 'none',
axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1,
size = 4),
axis.title = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
plot.background = element_blank(),
panel.grid = element_blank(),
plot.title = element_text(size = 5, margin=margin(0,0,0,0))) +
ggtitle(str_c(reg, '\nn=', pull(filter(regions.tally, region==reg))))
}
gg.accuracy <- region.plot.order %>%
set_names %>%
map(plot.accuracy, slide.regions, slide, modis)
cowplot::plot_grid(plotlist = gg.accuracy, nrow = 1)Climatology by region
plot.climatology <- function (reg, precip.monthly, slide.regions, med) {
precip.seasonal <- precip.monthly %>%
group_by(OBJECTID, year, season) %>%
summarize(precip.mm = sum(precip.mm))
ymin <- 1
ymax <- max(precip.seasonal$precip.mm, na.rm=T)
precip.seasonal %>%
dplyr::right_join(slide.regions %>% filter(region==reg)) %>%
drop_na(season) %>%
ggplot() +
geom_violin(aes(y = precip.mm, x = season, group = season),
fill='black', color=NA, adjust = 2) +
geom_hline(aes(yintercept = med), color = 'grey', size = .3)+
scale_y_log10(breaks = c(1, 10, 100, 1000)) +
scale_x_discrete(labels = c('Spring' = 'MAM', 'Summer' = 'JJA',
'Fall' = 'SON', 'Winter' = 'DJF'),
expand = ggplot2::expand_scale(add = c(1, 0))) +
coord_cartesian(expand=F, ylim = c(ymin, ymax)) +
theme_bw() +
theme(legend.position = 'none',
axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 4,
margin=margin(0,0,0,0)),
axis.text.y = element_text(size = 4, hjust = 1, vjust = 0.5),
axis.title = element_blank(),
plot.background = element_blank(),
panel.grid = element_blank(),
plot.title = element_text(size = 5, margin=margin(0,0,0,0))) +
ggtitle(str_c(reg, '\nn=', pull(filter(regions.tally, region==reg))))
}
gg.climatology <- region.plot.order %>%
purrr::set_names() %>%
purrr::map(plot.climatology,
precip.chirps.ww.monthly,
slide.regions,
precip.chirps.ww.seasonal.median)
cowplot::plot_grid(plotlist = gg.climatology, nrow = 1)Full map figure
coord_america <- coord_equal(xlim = c(-140, -50), ylim = c(-8, 53), expand = F)
coord_asia <- coord_equal(xlim = c(50, 140), ylim = c(-12, 49), expand = F)
# Clip maps
gg.america.region <- gg.region + theme_panel + coord_america
gg.asia.region <- gg.region + theme_panel + coord_asia
gg.america.burned <- gg.burned + theme_panel + coord_america
gg.asia.burned <- gg.burned + theme_panel + coord_asia
gg.america.precip <- gg.precip + theme_panel + coord_america
gg.asia.precip <- gg.precip + theme_panel + coord_asia
# Pull legends
legend.regions <- cowplot::get_legend(gg.region)
legend.burned <- cowplot::get_legend(gg.burned)
legend.precip <- cowplot::get_legend(gg.precip)
# Region midpoints for arrows
region.centroids <- slide.regions %>%
left_join(slide) %>%
group_by(region) %>%
summarize(y = mean(latitude), x = mean(longitude))# Add insets
inset.w <- .2
inset.h <- .41
inset.arrow <- grid::arrow(length = unit(0.03, "npc"), type = 'open')
# Regions
gg.america.region.inset <- cowplot::ggdraw(gg.america.region) +
# Arrows
cowplot::draw_grob(grid::grid.lines(x = c(.1, .25), y = c(.8, .8),
arrow=inset.arrow)) +
cowplot::draw_grob(grid::grid.lines(x = c(.1, .275), y = c(.4, .65),
arrow=inset.arrow)) +
cowplot::draw_grob(grid::grid.lines(x = c(.4, .4), y = c(.375, .575),
arrow=inset.arrow)) +
cowplot::draw_grob(grid::grid.lines(x = c(.825, .675), y = c(.675, .3),
arrow=inset.arrow)) +
# Inset plots
cowplot::draw_plot(gg.accuracy$`Pacific Northwest`,
x = 0.01, y = 1,
width = inset.w, height = inset.h, vjust = 1) +
cowplot::draw_plot(gg.accuracy$California,
x = 0.01, y = .6,
width = inset.w, height = inset.h, vjust = 1) +
cowplot::draw_plot(gg.accuracy$`Intermountain West`,
x = .275, y = .4,
width = inset.w, height = inset.h, vjust = 1) +
cowplot::draw_plot(gg.accuracy$`Central America`,
x = .775, y = .8,
width = inset.w, height = inset.h, vjust = 1)
gg.asia.region.inset <- cowplot::ggdraw(gg.asia.region) +
#Arrows
cowplot::draw_grob(grid::grid.lines(x = c(.3, .4), y = c(.3, .6),
arrow=inset.arrow)) +
cowplot::draw_grob(grid::grid.lines(x = c(.825, .8), y = c(.55, .4),
arrow=inset.arrow)) +
# Inset plots
cowplot::draw_plot(gg.accuracy$Himalayas,
x = .15, y = 0.4,
width = inset.w, height= inset.h, vjust = 1) +
cowplot::draw_plot(gg.accuracy$`Southeast Asia`,
x = .775, y = .935,
width = inset.w, height = inset.h, vjust = 1)# Fire
gg.america.burned.inset <- cowplot::ggdraw(gg.america.burned) +
cowplot::draw_plot(gg.burned.fraction$`Pacific Northwest`,
x = 0, y = 1,
width = inset.w, height = inset.h, vjust = 1) +
cowplot::draw_plot(gg.burned.fraction$California,
x = 0, y = .6,
width = inset.w, height = inset.h, vjust = 1) +
cowplot::draw_plot(gg.burned.fraction$`Intermountain West`,
x = .275, y = .4,
width = inset.w, height = inset.h, vjust = 1) +
cowplot::draw_plot(gg.burned.fraction$`Central America`,
x = .775, y = .8,
width = inset.w, height = inset.h, vjust = 1)
gg.asia.burned.inset <- cowplot::ggdraw(gg.asia.burned) +
cowplot::draw_plot(gg.burned.fraction$Himalayas,
x = .15, y = 0.4,
width = inset.w, height= inset.h, vjust = 1) +
cowplot::draw_plot(gg.burned.fraction$`Southeast Asia`,
x = .775, y = .935,
width = inset.w, height = inset.h, vjust = 1)
# Precipitation
inset.w.p = .23
inset.h.p = .41
gg.america.precip.inset <- cowplot::ggdraw(gg.america.precip) +
cowplot::draw_plot(gg.climatology$`Pacific Northwest`,
x = 0, y = 1.01,
width = inset.w.p, height = inset.h.p, vjust = 1) +
cowplot::draw_plot(gg.climatology$California,
x = .085, y = .64,
width = inset.w.p, height = inset.h.p, vjust = 1) +
cowplot::draw_plot(gg.climatology$`Intermountain West`,
x = .245, y = .395,
width = inset.w.p, height = inset.h.p, vjust = 1) +
cowplot::draw_plot(gg.climatology$`Central America`,
x = .71, y = .8,
width = inset.w.p, height = inset.h.p, vjust = 1)
gg.asia.precip.inset <- cowplot::ggdraw(gg.asia.precip) +
cowplot::draw_plot(gg.climatology$Himalayas,
x = .15, y = 0.4,
width = inset.w.p, height= inset.h.p, vjust = 1) +
cowplot::draw_plot(gg.climatology$`Southeast Asia`,
x = .7, y = .935,
width = inset.w.p, height = inset.h.p, vjust = 1)panel.labels <- tibble(
label = as_factor(c('Location accuracy distribution\nfor burned and unburned sites by region',
'Locations of burned sites and\nregional distribution of burned fraction',
str_c('Landslide-triggering precipitation percentile\nand',
'regional seasonal precipitation (mm)')))
)
gg.panel.labels <- ggplot(panel.labels) +
facet_wrap(~label, ncol = 1, switch = T) +
theme(panel.background = element_blank(),
strip.text = element_text(face="bold", size=7))
gg.panel.maps <- cowplot::plot_grid(
cowplot::plot_grid(gg.america.region.inset, gg.asia.region.inset,
labels = c('a', 'b'), label_x = 1, hjust = 1, vjust = 1),
legend.regions,
cowplot::plot_grid(gg.america.burned.inset, gg.asia.burned.inset,
labels = c('c', 'd'), label_x = 1, hjust = 1, vjust = 1),
legend.burned,
cowplot::plot_grid(gg.america.precip.inset, gg.asia.precip.inset,
labels = c('e', 'f'), label_x = 1, hjust = 1, vjust = 1),
legend.precip,
ncol = 1, rel_heights = c(1, .3, 1, .2, 1, .3))
gg.panel <- cowplot::plot_grid(gg.panel.labels, gg.panel.maps,
nrow = 1, rel_widths = c(.5, 6.5))
gg.panel
ggsave(filename = "../figures/panel_map.png", width=7, height=8)Fire timing
This figure shows the timing of fires relative to landslides for all events that took place at a burned site.
calc.fire.timing <- function (modis.date, modis, slide) {
modis.date %>%
dplyr::filter(burn_days_before >= 0) %>%
dplyr::group_by(OBJECTID) %>%
dplyr::filter(burn_days_before == min(burn_days_before)) %>%
right_join(modis %>% dplyr::filter(burned=='burned'), by='OBJECTID') %>%
inner_join(slide.regions, by='OBJECTID') %>%
left_join(slide %>% dplyr::select(OBJECTID, event_yday, season), by='OBJECTID') %>%
dplyr::mutate(burn_yday = yday(burn_date),
burn_season = factor(burn_season,
levels = c('Spring', 'Summer', 'Fall', 'Winter'),
ordered = T))
}
get.data('fire.timing', overwrite = F,
calc.fire.timing,
modis.date, modis, slide)Most landslides take place within one year of the fire. California stands out as having over a quarter of the landslides taking place more than a year after the fire.
fire.timing %>%
dplyr::mutate(gt.1year = ifelse(burn_days_before > 2*365,
'gt.1year', 'lt.1year')) %>%
dplyr::group_by(region, gt.1year) %>%
dplyr::tally() %>%
tidyr::spread(gt.1year, n, drop=F) %>%
dplyr::mutate(prop.gt.1year = gt.1year / (gt.1year + lt.1year))fire.to.landslide.x.breaks <- seq(season.breaks[4] - 4 * 365, season.breaks[3],
365/4)
fire.to.landslide.x.labels <- map(c(str_c('-', c(3, 2, 1), 'y '), ''),
~c(str_c(., 'Winter'),
'Spring', 'Summer', 'Fall')) %>% as_vector
fire.facet.labels <- fire.timing %>%
dplyr::group_by(region) %>%
dplyr::mutate(wrap_event_yday = ifelse(event_yday > season.breaks['winter'],
event_yday - 365, event_yday)) %>%
dplyr::ungroup() %>%
dplyr::mutate(ymin = min(wrap_event_yday - burn_days_before),
ymax = max(wrap_event_yday)) %>%
dplyr::group_by(region) %>%
dplyr::summarize(y = n() * 0.8, x = min(ymin) + (max(ymax) - min(ymin)) * 0.1) %>%
dplyr::mutate(label = c('a', 'b', 'c', 'd', 'e', 'f', 'g'))
fire.timing %>%
dplyr::group_by(region) %>%
dplyr::mutate(wrap_event_yday = ifelse(event_yday > season.breaks['winter'],
event_yday - 365, event_yday)) %>%
dplyr::arrange(wrap_event_yday - burn_days_before) %>%
dplyr::mutate(i = row_number()) %>%
ggplot() +
facet_wrap(vars(region), scales = 'free_y', ncol = 2, shrink=F) +
geom_segment(aes(x = wrap_event_yday - burn_days_before, xend = wrap_event_yday,
y = i, yend = i,
color = burn_season)) +
geom_point(aes(x = wrap_event_yday, y = i, shape='Landslide'), size=0.5) +
geom_vline(data = tibble(x=seq(334-365*3, 365, 365)),
aes(xintercept = x, size = 'start_of_winter')) +
geom_vline(aes(xintercept = 334 - 365, size = 'year_of')) +
geom_rug(aes(y = i, color = burn_season)) +
geom_rug(aes(x = event_yday - burn_days_before), sides = 't', length = unit(0.05, "npc")) +
geom_label(data = fire.facet.labels,
aes(x = x, y = y, label = label)) +
theme_bw() +
scale_x_continuous('Landslides and preceding fires relative to the year of the landslide',
minor_breaks = NULL,
breaks = as_vector(fire.to.landslide.x.breaks),
labels = fire.to.landslide.x.labels) +
scale_color_manual('Fire season:',
values = c('Spring'='#19E667', 'Summer'='#E61998',
'Fall'='#E6CD19', 'Winter'='#1932E6'),
guide = guide_legend(nrow=1, order = 2)) +
scale_y_continuous('Post-wildfire landslides ordered by fire timing relative to landslide',
expand = c(0.1, 0.1)) +
scale_size_manual('Relative timing:',
values = c('year_of' = 1, 'start_of_winter' = 0.5),
labels = c('year_of' = 'Start of year of landslide',
'start_of_winter' = 'Start of winter (December 1)'),
guide = guide_legend(nrow=2, order = 3)) +
guides(shape = guide_legend(title = NULL, override.aes = c(size = 1.5), order = 1)) +
theme(panel.grid.major.x = element_line(color = 'grey', size = 0.5),
axis.text.x = element_text(hjust = 1, vjust = 1, angle = 90),
axis.ticks = element_blank(),
axis.text.y = element_blank(),
legend.position = c(.75, .04),
legend.direction = 'vertical',
legend.box = 'vertical',
legend.margin = unit(5, 'pt'),
legend.box.margin = margin(5, 5, 5, 5),
legend.box.background = element_rect(color = 'black', size = 2),
legend.key.height = unit(10, 'pt')
)
ggsave('../figures/fire_to_landslide.png', width = 6.5, height = 8)Analysis
CHIRPS vs Daymet
To validate the exclusion of zero-precipitation landslides, a comparison of CHIRPS and Daymet in California and Nevada only. Many of the zero-precipitation landslides measured by CHIRPS are not the same as the ones measured by Daymet, suggesting that these are measurement errors.
Approximately 14% of the landslides in California and Nevada were excluded on this basis.
join.chirps.daymet <- function (precip.chirps, precip.daymet) {
precip.chirps %>%
dplyr::inner_join(precip.daymet %>% dplyr::select(-days_before, -yday),
by=c('OBJECTID', 'date', 'window'),
suffix=c('.chirps', '.daymet'))
}
get.data('precip.chirps.daymet.ca.event', join.chirps.daymet, overwrite = F,
precip.chirps.ca.event,
precip.daymet.ca.all)
get.data('precip.chirps.daymet.ca.all', join.chirps.daymet, overwrite = F,
precip.chirps.ca.all, precip.daymet.ca.all)
precip.chirps.ca.event %>%
filter(days_before==-1, window==7, !is.na(pctl)) %>%
group_by(precip = ifelse(pctl > 0, 'non.zero.precip', 'zero.precip')) %>%
tally() %>%
spread(key = precip, value = n) %>%
mutate(total = non.zero.precip + zero.precip,
percent.zero.precip = zero.precip / total)precip.chirps.daymet.ca.all %>%
tidyr::drop_na(pctl.chirps, pctl.daymet) %>%
ggplot() +
geom_point(aes(x=pctl.chirps, y=pctl.daymet, color = 'all'),
alpha=0.2,shape = 16) +
geom_point(data = precip.chirps.daymet.ca.event %>%
dplyr::filter(window == 7, days_before == -1 ) %>%
dplyr::mutate(to.remove = ifelse(pctl.chirps==0, 'exclude', 'keep')),
aes(x=pctl.chirps, y=pctl.daymet, color=to.remove)) +
coord_equal() +
theme_bw() +
scale_color_manual('',
values = c('keep'='black', 'exclude'='blue', 'all'='gray'),
labels = c('keep'='Included\nday-of-landslide\npreciptiation',
'exclude'='Excluded\nday-of-landlside\nprecipitation',
'all'='Precipitation\nfrom all dates')) +
guides(color = guide_legend(keyheight = grid::unit(4, 'char'))) +
labs(x='CHIRPS 7-day Precipitation Percentile',
y='Daymet 7-day Precipitation Percentile')
ggsave('../figures/chirps_vs_daymet.png', width=6, height=4, dpi=600)Precipitation percentile wilcox test
Comparison of antecedent precipitation time series among burned and unburned sites by region. Mann-Whitney tests are used to establish whether or not differences are statistically significant.
calc.burn.wilcox <- function (precip, pre=30, post=10, alternative='two.sided',
source=NA, groupby=c('burned', 'days_before')) {
precip %>%
dplyr::filter(days_before <= pre, days_before >= -post) %>%
dplyr::group_by_at(groupby) %>%
calc.wilcox(alternative, source)
}
calc.wilcox <- function (precip, alternative='two.sided', source=NA) {
wilcox <- precip %>%
tidyr::nest() %>%
spread(key = 'burned', value = 'data') %>%
dplyr::mutate(wilcox = map2(burned, unburned,
~wilcox.test(.x$pctl, .y$pctl,
alternative=alternative))) %>%
dplyr::mutate(wilcox.glance = map(wilcox, glance)) %>%
tidyr::unnest(wilcox.glance) %>%
dplyr::mutate(signif.95 = p.value <= 0.05) %>%
gather(key = 'burned', value = 'data', burned, unburned) %>%
tidyr::unnest(data)
if (!is.na(source)) wilcox <- wilcox %>% dplyr::mutate(label=source)
wilcox
}
magnitude.regional.data <- function (slide.regions, precip,
pre = 30, post = 10, w = 7) {
precip %>%
dplyr::filter(days_before < pre, days_before > -post, window==w) %>%
right_join(slide.regions, by = 'OBJECTID') %>%
dplyr::left_join(modis, by = 'OBJECTID') %>%
dplyr::filter(pctl > 0) %>%
droplevels() %>%
dplyr::select(region, burned, days_before, pctl) %>%
calc.burn.wilcox(pre = pre, post = post, alternative = 'less',
groupby = c('region', 'burned', 'days_before'))
}
get.data('burn.wilcox.regions.event', overwrite = F,
magnitude.regional.data,
slide.regions,
precip.chirps.ww.event)facet.levels <- c('Rolling cumulative\nprecipitation percentile',
'Mann-Whitney p-value')
plot.burn.wilcox <- function (burn.wilcox.regions.event, slide.regions) {
region.tally <- slide.regions %>%
left_join(modis) %>%
dplyr::group_by(burned=burn.cumulative > 0, region) %>%
tally
burn.wilcox.regions.event %>%
dplyr::mutate(signif.95 = ifelse(signif.95, 'significant', 'not.significant')) %>%
ggplot() +
facet_wrap(region~., ncol=1, strip.position = 'right') +
geom_boxplot(aes(x = days_before,
group=interaction(days_before, burned),
y = pctl,
fill = burned,
alpha = interaction(signif.95, burned)),
outlier.shape=21) +
annotate("rect", xmin = -0.5, xmax = 0.5, ymin = 0, ymax = 1, alpha = .4) +
geom_label(data=region.tally,
aes(x=-15, y=0.1+0.25*burned,
label=str_c(ifelse(burned, 'Burned: ', 'Unburned: '), n)),
hjust='right', label.size = 0) +
geom_label(data = tibble(region = region.plot.order,
label = c('a', 'b', 'c', 'd', 'e', 'f', 'g')) %>%
dplyr::mutate(region = factor(region, region.plot.order, ordered = T)),
aes(x = -15, y = 0.85, label = label), hjust = 'right') +
scale_x_reverse(limits = c(31, -16), expand = expansion()) +
scale_fill_manual('',
values = c('burned' = '#d95f02', unburned = '#7570b3'),
labels = c('burned'='Burned', 'unburned'='Unburned'),
guide = guide_legend(direction = 'horizontal')) +
scale_alpha_manual('Statistical significance at\n95% confidence',
values=c('not.significant.burned' = 0.4,
'not.significant.unburned' = 0.4,
'significant.burned' = 1,
'significant.unburned' = 1),
labels=c('not.significant.burned' = '',
'not.significant.unburned' = 'Not Significant',
'significant.burned' = '',
'significant.unburned' = 'Significant'),
breaks=c('not.significant.burned',
'not.significant.unburned',
'significant.burned',
'significant.unburned'),
guide=guide_legend(override.aes=list(fill=c('not.significant.burned' = '#d95f02',
'not.significant.unburned' = '#7570b3',
'significant.burned' = '#d95f02',
'significant.unburned' = '#7570b3')),
direction = 'horizontal')) +
theme_bw() +
labs(x='Days before the event', y='Rolling cumulative precipitation percentile') +
theme(legend.position = 'bottom', legend.direction = 'vertical')
}
burn.wilcox.regions.event %>%
plot.burn.wilcox(slide.regions)
ggsave('../figures/magnitude_regional.png', width = 11, height = 11)Precipitation (mm) bootstrap sampling
Bootstrap samples are used here to develop baseline precipitation values against which to compare landslide-triggering precipitation. This method is an improvement over the original Wilcox tests because sample-size disparities between regions are reduced.
Bootstrap sampling
sample.bootstrap <- function(precip.subset, precip.all, n, samples, burned, region) {
precip.all %>%
dplyr::filter(burned==burned, region==region) %>%
dplyr::mutate(precip.yday = yday(date)) %>%
# dplyr::select precipitation from same location and date of year as subset
dplyr::inner_join(precip.subset %>%
dplyr::select(OBJECTID, event_yday),
by=c('OBJECTID'='OBJECTID',
'precip.yday'='event_yday')) %>%
# Sample precipitation matching the site/day in bulk
dplyr::group_by(OBJECTID) %>%
tidyr::drop_na(pctl) %>%
dplyr::sample_n(n, replace = T) %>%
# Randomly split into individual samples
dplyr::mutate(i = sample(1:samples, n, replace = T ))
}
calc.burn.bootstrap <- function (precip.event, precip.all,
samples=100, per.subset=10000,
pre=30, post=10, w=7) {
groupby = c('region', 'burned', 'days_before')
precip.all <- precip.all %>%
dplyr::filter(window==w) %>%
dplyr::select(-window)
precip.event %>%
rename(precip.date = date, precip.yday = yday) %>%
dplyr::filter(days_before <= pre, days_before >= -post, window==w) %>%
dplyr::select_at(c(groupby, 'OBJECTID', 'event_yday')) %>%
# Tally each subset/site/day combination
dplyr::group_by_at(groupby) %>%
add_tally(name = 'sites.in.subset') %>%
# Determine number to draws from each site that will
# roughly match sample sizes among subsets
# Round up - number of draws must be an integer
dplyr::mutate(n.per.site = ceiling(per.subset / sites.in.subset)) %>%
# Take bootstrap samples for each subset/site/day combination
dplyr::group_by_at(c(groupby, 'n.per.site', 'sites.in.subset')) %>%
tidyr::nest() %>%
dplyr::mutate(bootstrap.samples = pmap(list(data, n.per.site,
burned, region),
~sample.bootstrap(..1, precip.all,
..2, samples,
..3, ..4)))
}join.details <- function (precip, slide, slide.regions, modis) {
precip %>%
dplyr::inner_join(slide %>% dplyr::select(OBJECTID, event_yday),
by='OBJECTID') %>%
dplyr::right_join(slide.regions, by = 'OBJECTID') %>%
dplyr::left_join(modis %>% dplyr::select(OBJECTID, burned),
by='OBJECTID')
}
get.bootstrap <- function (reg='All landslides',
name = 'precip.chirps.bootstrap.all',
ovw=F) {
get.data(name, overwrite = ovw,
calc.burn.bootstrap,
precip.chirps.ww.event %>%
join.details(slide, slide.regions %>% filter(region == reg), modis),
precip.chirps.ww.all %>%
join.details(slide, slide.regions %>% filter(region == reg), modis))
}
get.bootstrap('California', 'precip.chirps.bootstrap.cal', ovw = F)
get.bootstrap('Intermountain West', 'precip.chirps.bootstrap.imw', ovw = F)
get.bootstrap('Pacific Northwest', 'precip.chirps.bootstrap.pnw', ovw = F)
get.bootstrap('Central America', 'precip.chirps.bootstrap.cam', ovw = F)
get.bootstrap('Himalayas', 'precip.chirps.bootstrap.him', ovw = F)
get.bootstrap('Southeast Asia', 'precip.chirps.bootstrap.sea', ovw = F)
get.bootstrap(ovw = F)
precip.chirps.bootstrap <- precip.chirps.bootstrap.all %>%
dplyr::bind_rows(precip.chirps.bootstrap.cal) %>%
dplyr::bind_rows(precip.chirps.bootstrap.imw) %>%
dplyr::bind_rows(precip.chirps.bootstrap.pnw) %>%
dplyr::bind_rows(precip.chirps.bootstrap.cam) %>%
dplyr::bind_rows(precip.chirps.bootstrap.him) %>%
dplyr::bind_rows(precip.chirps.bootstrap.sea)
precip.chirps.bootstrap %>%
dplyr::filter(days_before==0) %>%
dplyr::select(-data) %>%
dplyr::mutate(site.total = map_int(bootstrap.samples, ~nrow(ungroup(.))),
n.per.site.real = site.total / sites.in.subset,
avg.sample.size = map_dbl(bootstrap.samples,
~group_by(., i) %>%
summarize(n = n()) %>%
pull(n) %>% mean)) %>%
dplyr::select(-bootstrap.samples)Summarize bootstrap samples
mann.whitney.p.value <- function (subset, bootstrap) {
bootstrap %>%
dplyr::group_by(i) %>%
tidyr::nest() %>%
# test if bootstrap samples are less than pre-landslide precipitation
dplyr::mutate(wilcox = map(data,
~wilcox.test(pull(., pctl),
pull(subset, pctl),
alternative='less'))) %>%
dplyr::mutate(wilcox.glance = map(wilcox, glance)) %>%
dplyr::select(-data, -wilcox) %>%
tidyr::unnest(wilcox.glance)
}subset.precip <- function (precip, regions, modis, w=7, pre=30, post=10) {
precip %>%
dplyr::filter(window==7, days_before <= pre, days_before >= -post) %>%
dplyr::select(OBJECTID, pctl, days_before) %>%
dplyr::right_join(regions, by='OBJECTID') %>%
dplyr::inner_join(modis, by='OBJECTID') %>%
dplyr::group_by(region, burned, days_before) %>%
tidyr::nest()
}
get_subset <- function (reg, bur, day, subsets) {
row <- subsets %>%
dplyr::filter(region==as.character(reg),
burned==bur,
days_before==day)
row$data[[1]]
}
get.data('precip.chirps.ww.event.subsets',
subset.precip,
precip.chirps.ww.event,
slide.regions,
modis)summarize.bootstrap <- function (bootstrap.df, subsets, objective.function) {
bootstrap.df %>%
dplyr::mutate(objective = purrr::pmap(list(region,
burned,
days_before,
bootstrap.samples),
~objective.function(get_subset(..1,
..2,
..3,
subsets),
..4))) %>%
dplyr::select(-bootstrap.samples) %>%
tidyr::unnest(objective) %>%
dplyr::ungroup()
}get.data('precip.chirps.bootstrap.wilcox', overwrite = F,
summarize.bootstrap,
precip.chirps.bootstrap %>% dplyr::select(-data),
precip.chirps.ww.event.subsets,
mann.whitney.p.value)
precip.chirps.bootstrap.wilcox %>%
dplyr::select(region, burned, days_before, p.value)The bootstrap plot also shows the day-of-landslide kernel density estimates for all bootstrap samples against the landslide-triggering precipitation to illustrate the method. ### Plot bootstrap sampling
plot.burn.bootstrap <- function (bootstrap, var, name) {
bootstrap %>%
ggplot() +
facet_wrap(.~region, ncol=1, strip.position='right') +
geom_vline(aes(xintercept=0)) +
geom_boxplot(aes(x = days_before,
group = interaction(days_before, burned),
y = {{ var }},
fill = burned,
color = burned),
outlier.shape=21,
alpha = 0.6,
position = position_dodge(width=0.4, preserve = 'total')) +
annotate("rect", xmin = -11, xmax = -0.1, ymin = 0.05, ymax = .9999,
fill = 'white', alpha = 0.4) +
annotate("rect", xmin = 0.1, xmax = 31, ymin = 0.05, ymax = .9999,
fill = 'white', alpha = 0.4) +
geom_hline(aes(yintercept=0.05), linetype='dashed') +
geom_label(data = tibble(region = region.plot.order,
label = c('a', 'b', 'c', 'd', 'e', 'f', 'g')) %>%
dplyr::mutate(region = factor(region, region.plot.order, ordered = T)),
aes(x = -11.5, y = 0.9999, label = label),
hjust = 'right', vjust = 'bottom') +
scale_x_reverse(expand = expansion(), limits = c(31, -12)) +
scale_y_continuous(trans = probit_trans(),
breaks = bootstrap.y.breaks,
labels = bootstrap.y.labels) +
scale_fill_burned +
scale_color_burned +
coord_trans(y = 'reverse',
ylim = c(10^-35, .9999)) +
theme_bw() +
labs(x='Days before the landslide', y=name) +
theme(legend.position = 'bottom',
panel.grid.minor.y = element_blank())
}
bootstrap.y.breaks <- c(10^-30, 10^-10, 0.05, 0.5, 0.9999)
bootstrap.y.labels <- c('1e-30', '1e-10', '0.05', '0.5', '1')
plot.bootstrap.wilcox <- precip.chirps.bootstrap.wilcox %>%
plot.burn.bootstrap(p.value,
'p-value of Wilcox tests between bootstrap samples and pre-landslide precipitation')
plot.bootstrap.wilcoxplot.bootstrap.sample.kde <- precip.chirps.bootstrap %>%
dplyr::filter(days_before==0) %>%
dplyr::select(-data) %>%
tidyr::unnest() %>%
ggplot() +
facet_grid(cols = vars(burned), rows = vars(region), scales = 'free_y') +
geom_density(aes(x = pctl, group = i, color = burned, fill = burned),
alpha = 0.01) +
geom_density(data = precip.chirps.ww.event %>%
filter(days_before==0) %>%
right_join(slide.regions, by = 'OBJECTID')%>%
inner_join(modis, by='OBJECTID'),
aes(x = pctl, color = 'day_of', fill = 'day_of'), size = 1) +
geom_label(data = cross_df(list(region = region.plot.order,
burned = c('burned', 'unburned'))) %>%
dplyr::mutate(label = c('h', 'i', 'j', 'k', 'l', 'm', 'n',
'o', 'p', 'q', 'r', 's', 't', 'u')) %>%
dplyr::mutate(region = factor(region, region.plot.order, ordered = T)),
aes(x = .95, y = .05, label = label),
hjust = 'right', vjust = 'bottom') +
scale_fill_manual('',
breaks = c('day_of', 'burned', 'unburned'),
values = c('day_of' = NA, burned.colors),
labels = c('day_of' = 'Day of landslide precipitation',
'burned' = 'Bootstrap samples, burned sites',
'unburned' = 'Bootstrap samples, unburned sites'),
guide = guide_legend(override.aes = list(alpha = 0.5))) +
scale_color_manual('',
breaks = c('day_of', 'burned', 'unburned'),
values = c('day_of' = 'black', burned.colors),
labels = c('day_of' = 'Day of landslide precipitation',
'burned' = 'Bootstrap samples, burned sites',
'unburned' = 'Bootstrap samples, unburned sites'),
guide = guide_legend(override.aes = list(alpha = 0.5))) +
scale_x_continuous(breaks = seq(0, 1, 0.25),
labels = c('0', '0.25', '0.5', '0.75', '1')) +
labs(y = 'Kernel density estimate', x = 'Precipitation percentile') +
theme_bw() +
theme(legend.position = 'bottom',
legend.spacing.x = unit(15, "pt"),
axis.text.y = element_blank(),
axis.text.x = element_text(),
axis.ticks.y = element_blank(),
strip.text.x = element_blank(),
panel.spacing.x = unit(15, 'pt')) +
coord_cartesian(expand = F)
plot.bootstrap.sample.kdeplot.bootstrap.sample.kde.legend <- get_legend(plot.bootstrap.sample.kde)
plot.bootstrap.wilcox.legend <- get_legend(plot.bootstrap.wilcox)
cowplot::plot_grid(cowplot::plot_grid(plot.bootstrap.wilcox +
theme(strip.text = element_blank(),
plot.margin = unit(c(5,5,5,5), 'pt'),
legend.position = 'none',
legend.direction = 'horizontal'),
plot.bootstrap.sample.kde +
theme(plot.margin = unit(c(5,5,5,2), 'pt'),
legend.position = 'none'),
rel_widths = c(11, 4)),
cowplot::plot_grid(plot.bootstrap.wilcox.legend,
plot.bootstrap.sample.kde.legend,
rel_widths = c(3, 7)),
nrow = 2, rel_heights = c(10, 0.5))
ggsave('../figures/precip_bootstrap_wilcox+kde.png', width=13, height=11)Seasonal density and precipitation frequency
Seasonal density
The kernel density estimate of the annual distribution of landslides is used to show shifts in seasonality relative to the day of the year.
slide.season <- slide %>%
dplyr::select(OBJECTID, location_accuracy, season, event_yday) %>%
dplyr::right_join(slide.regions, by='OBJECTID') %>%
dplyr::left_join(modis, by='OBJECTID')
plot.seasonality <- bind_rows(slide.season %>%
dplyr::mutate(event_yday = event_yday - 365),
slide.season,
slide.season %>%
dplyr::mutate(event_yday = event_yday + 365)) %>%
dplyr::group_by(region, burned) %>%
tidyr::drop_na(burned) %>%
tidyr::nest() %>%
dplyr::mutate(kde = map(data,
function (df) density(df$event_yday,
from=-100, to=450,
n = 512, bw = 21)),
kde = map(kde, function (dens) tibble(event_yday = dens$x,
density = dens$y))) %>%
tidyr::unnest(kde) %>%
dplyr::filter(event_yday >= 0, event_yday < 365) %>%
ggplot() +
facet_wrap(region~., ncol = 1, strip.position = 'right') +
geom_area(aes(x = event_yday, y = density, group = burned, fill=burned),
trim = T, alpha = 0.6, color = 'black',
position = position_identity()) +
geom_text(data = tibble(text = c('h', 'i', 'j', 'k', 'l', 'm', 'n'),
hjust = c(-3.5, -10, -10, -4, -10, -2, -3.5),
region = region.plot.order,
x = -2.1*365, y = 0.50) %>%
mutate(region = factor(region, region.plot.order, ordered = T)),
aes(x = 0, y = .0035, label = text, hjust = hjust), size = 7) +
coord_polar() +
scale_fill_burned +
scale_x_continuous(breaks = (season.breaks + 365/8) %% 365,
labels = c('W', 'Sp', 'Su', 'F'),
minor_breaks = season.breaks) +
labs(y = '', x = ' ', fill = 'Season') +
theme_bw() +
theme(axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
axis.text.x = element_text(angle = 0, vjust = 1),
legend.position = 'bottom',
panel.grid.minor.x = element_line(color = 'black'))
plot.seasonalityPrecipitation frequency
Smoothed precipitation frequency is used to show shifts in the annual distribution of landslides relative to the annual precipitation seasonality.
n.years <- 2
calc.prop <- function (col) sum(!(col %in% 0)) / n()
calc.yday.frequency <- function (precip, slide, slide.regions, modis) {
precip %>%
dplyr::mutate(precip_yday = yday(date)) %>%
dplyr::select(OBJECTID, precip_yday, pctl) %>%
dplyr::inner_join(slide %>% dplyr::select(OBJECTID, event_yday),
by='OBJECTID') %>%
dplyr::mutate(ydays_before = event_yday - precip_yday) %>%
dplyr::select(OBJECTID, ydays_before, pctl) %>%
dplyr::mutate(ydays_before = ydays_before %% 365) %>%
dplyr::right_join(slide.regions, by = 'OBJECTID') %>%
dplyr::left_join(modis %>% dplyr::select(OBJECTID, burned),
by='OBJECTID') %>%
dplyr::group_by(region, burned, ydays_before) %>%
dplyr::summarize_at(vars(pctl), list(precip.yday.freq = calc.prop))
}
get.data('precip.chirps.yday.frequency',
calc.yday.frequency,
precip.chirps.ww.all,
slide,
slide.regions,
modis)
calc.frequency <- function (precip, precip.yday.freq,
slide, slide.regions, modis) {
precip.chirps.ww.all %>%
dplyr::transmute(OBJECTID, date = lubridate::as_date(date), pctl) %>%
dplyr::inner_join(slide %>% dplyr::select(OBJECTID, event_date),
by='OBJECTID') %>%
dplyr::mutate(days_before = as.integer(event_date - date)) %>%
dplyr::select(OBJECTID, days_before, pctl) %>%
dplyr::filter(days_before < n.years * 365,
days_before > -n.years * 365) %>%x
dplyr::right_join(slide.regions, by = 'OBJECTID') %>%
dplyr::left_join(modis %>% dplyr::select(OBJECTID, burned),
by='OBJECTID') %>%
dplyr::group_by(region, burned, days_before) %>%
dplyr::summarize_at(vars(pctl), list(precip.day.freq = calc.prop)) %>%
dplyr::mutate(ydays_before = days_before %% 365) %>%
dplyr::left_join(precip.yday.freq, by=c('region', 'ydays_before', 'burned')) %>%
dplyr::mutate(precip.freq.norm = precip.day.freq / precip.yday.freq,
precip.freq.anomaly = precip.day.freq - precip.yday.freq,
precip.freq.lt.anomaly = precip.day.freq - mean(precip.day.freq)) %>%
dplyr::ungroup()
}
get.data('precip.chirps.frequency', overwrite = F,
calc.frequency,
precip.chirps.ww.all,
precip.chirps.yday.frequency,
slide,
slide.regions,
modis)plot.frequency <- function (data, column, y.lable, n.years, groupby=c('burned')) {
data %>%
dplyr::group_by_at(groupby) %>%
dplyr::mutate(smooth = rollmean(!!as.name(column), 90, na.pad=TRUE)) %>%
ggplot() +
geom_line(aes(x = days_before, y = !!as.name(column),
group = burned, color = burned, size = 'daily')) +
geom_line(aes(x = days_before, y = smooth, group = burned, color = burned,
size = 'rolling')) +
geom_hline(data = data %>%
dplyr::group_by_at(groupby) %>%
summarize(mean = mean(!!as.name(column))),
aes(yintercept = mean, color=burned), linetype = 'dashed') +
geom_text(data = tibble(text = c('a', 'b', 'c', 'd', 'e', 'f', 'g'),
region = region.plot.order,
x = -1.95*365, y = 0.50),
aes(x = x, y = y, label = text), size = 7,
hjust = 'right') +
scale_size_manual('',
values = c('daily' = 0.1, 'rolling' = 0.75),
labels = c('daily' = 'Daily frequency',
'rolling' = '90-day rolling average')) +
scale_color_manual('',
values = c('burned' = '#d95f02', unburned = '#7570b3'),
labels = c('burned'='Burned', 'unburned'='Unburned'),
guide = guide_legend(override.aes = list(size = 2))) +
scale_x_reverse('Years before the landslide',
breaks = seq(-365 * n.years, 365 * n.years, 365),
labels = function (breaks) breaks / 365,
expand = expansion()) +
labs(y = y.lable) +
theme_bw() +
theme(legend.position = 'bottom')
}
plot.frequency.region <- precip.chirps.frequency %>%
plot.frequency('precip.freq.lt.anomaly',
'Precipitation frequency anomaly - long-term mean',
n.years=n.years,
groupby=c('burned', 'region')) +
facet_wrap(region~., ncol=1, strip.position = 'right')
plot.frequency.regioncowplot::plot_grid(
cowplot::plot_grid(
plot.frequency.region +
theme(strip.text.y = element_blank(),
legend.position = 'none',
plot.margin = unit(c(5, 1, 5, 5), 'pt')),
plot.seasonality +
theme(legend.position = 'none',
plot.margin = unit(c(5, 5, 12, 0), 'pt')),
rel_widths = c(10, 2.5)
),
get_legend(plot.frequency.region),
nrow = 2, rel_heights = c(10, 0.5)
)
ggsave('../figures/precip_frequency_seasonality.png', width=11, height=12)